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We investigate the effective conductivity (a e ) of a class of amorphous media defined by the level- 
cut of a Gaussian random field. The three point solid-solid correlation function is derived and utilised 
in the evaluation of the Beran-Milton bounds. Simulations are used to calculate cr e for a variety of 
fields and volume fractions at several different conductivity contrasts. Relatively large differences 
in <r e are observed between the Gaussian media and the identical overlapping sphere model used 
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I. INTRODUCTION 



The calculation of the effective transport properties of 
random composite media is important in many scientific 
and engineering applications pj . Several techniques (ef- 
fective medium approximations and cluster expansions) 
have been developed for predicting the effective proper- 
ties of such materials (briefly reviewed in ref. g]). How- 
ever difficulties encountered in such methods have pro- 
vided the impetus for the development of rigorous bounds 
|^-|| . Such bounds rely on statistical descriptions of the 
microstructure of the material which are available for rel- 
atively few classes of media. The advancement of com- 
puting technology f|-[lj] has also made direct simulation 
of effective properties feasible. It is the latter two ap- 
proaches which we shall discuss here, in the context of 
the effective conductivity of a three dimensional amor- 
phous isotropic two phase material. 

There has been significant advances in the evalua- 
tion of the Beran-Milton (BM) bounds in the past 
decade. The key parameter £i (or £ 2 = 1 — Ci) which 
incorporates microstructural information regarding the 
composite has been evaluated primarily for materials 
comprised of statistically independent cells jl3 ,fl4| o r dis- 
persions of regularly shaped inclusions |il^]llf 2yl7|] . The 
simulation of the effective conductivity of continuum ran- 
dom media is a computationally intensive process and has 
only recently been studied for the second class of materi- 
als ficjjllj . Such an approach provides a basis for testing 
the bounding theories and for generating outright predic- 
tions of the effective properties of composites. 

A class of materials which is not, in general, well de- 
scribed by cellular or particulate models is that of amor- 
phous composites. Such materials arise in certain alloys 
p8|,^9[ , microemulsions ]20| , ^l[ | and other systems p2[ . 
The model which best captures some of the salient fea- 



tures of such composites is the spatially uncorrelated pen- 
etrable sphere (or the identical overlapping sphere) model 
p3| . Due to the simplicity in evaluating the statistical 
correlation functions of such a material it has served as 
a useful 'model' amorphous medium |^J,|l5|,|5). However 
specific features of this model restrict its generality. The 
inclusion (sphere) phase and the matrix phase are topo- 
logically very different, the small scale structure of the 
phase boundaries is spherical and there are no long range 
correlations in the model. An alternative approach is 
to empirically measure the specific correlation functions 
of a sample and to apply the results in the evaluation 
of bounds (2^jl^,^,^|. This approach is complicated 
and subject to error. It is therefore interesting to seek a 
more complex model of amorphous composites, yet sim- 
ple enough so that the correlation functions can be cal- 
culated. 

Another method of modeling random composites is to 
define the interface between the phases as a level cut of 
some random fiel d p9| - p3[ | (see ref. (34|] for a review) . Re- 
cent progress |35| |37|] in the theory of interfaces of level- 
cut Gaussian random fields has made it possible to cal- 
culate the statistical information necessary for the eval- 
uation of the BM bounds. There is evidence that the 
Gaussian random interface model is a goo d approxima- 
tion to certain oil- water microemulsions fl38| , |39l , pi| and we 
conjecture that it is a reasonable model for amorphous 
alloys. 

In this paper we investigate the effective conductiv- 
ity of such media using the above mentioned bounding 
techniques and computer simulations. The results are 
compared with the previously studied models to demon- 
strate the differences that arise. The paper is organised 
as follows. In section [n] we describe the equations gov- 
erning the electric field in a composite medium and the 
bounds on the effective conductivity. In sections III & 
IV we derive the statistical correlation functions for the 
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random media and apply them in the calc ulation of the 
microstructure parameter. Sections |v|- |VI]| are concerned 
with the generation of the random materials, the simula- 
tion of the effective conductivity and comparison of the 
data with the bounds. 



II. BOUNDS ON THE EFFECTIVE PROPERTIES 
OF COMPOSITE MATERIALS 

The relationship between the current density j and the 
electric field E — — V0 is given by Ohms Law, 

j = -aV0. (1) 

Where, due to charge conservation, satisfies, 



X7 z (h = 



(2) 



throughout the material. At the boundary of different 
regions of the material with conductivities o\ and 02 we 
have, 



^1 = 4>2, criV^i.n = tT 2 V02-n 



(3) 



The effective conductivity is defined by a macroscopic 
form of Ohm's law, 



< erV0 > 
< V4> > ' 



(4) 



Now consider a composite material made up of two com- 
ponents with conductivities a% and 172 with volume frac- 
tions p and q = 1— p. The effective conductivity will then 
depend on the Cj, their respective volume fractions, and 
the spatial distribution (microstructure) of each phase 

The first bounds on a e were calculated by Wiener ]4l] ] 
who proved that < cr _1 > _1 < a e << a >. These bounds 
assume no details about the microstructure and are hence 
valid for a general composite. As more statistical infor- 
mation regarding the composite is included in the cal- 
culation of the bounds they become more restrictive. If 
the sample is assumed to be isotropic and macroscopi- 
cally homogeneous then the 2nd order bounds of Hashin 
and Sthrikman [S| are applicable. To distinguish between 
such materials the third order bounds of Beran are nec- 
essary. (The term nth order bounds refers to the fact 
that the bounds are exact to 0{u\ — 02)™). The Be- 
ran Q bounds were derived using variational principles 
and were subsequently simplified by Milton [Q. Follow- 
ing the notation of Milton we define < a >= pa\ + qa2, 
< a >= qai + pa,2 (interchanging p and q) and < a >(— 
Ci fl i + C2 a 2- Here a,; = <Ji or 1/cr.;. In these terms the 
lower bound on a e is, 
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< a 



> 



2pq(a 1 -a 2 ) 



2 < a- 1 > 
while the upper bound is, 



< cr 



(5) 



< a > -- 



pq(p x - <7 2 ) 



< a > +2 < cr >/ 



(6) 



The so called microstructure parameter £1 is given by a 
number of equivalent integrals , of which the for- 

mulation due to Brown Eot] is the best for our purposes, 



9 f^dr f^ds f 1 

2^ V T / duP2{u) x 
z pi Jo r Jo s J-i 



P3(r,s,t) 



P2{r)p2(s) 



(7) 



where t 2 = 



2rsu and P 2 (u) = (3u 2 - l)/2 is the 



Legendre polynomial of order 2. The functions p n are 
n-point solid-solid correlation functions (see section III) 
where the 'solid' is phase 1 and the 'void' is phase 2. 

As Milton notes these bounds converge when £i =0,1 
and are equal to one of the 2nd order Hashin-Sthrikman 
bounds in each case. An improved lower bound has been 
derived by Milton [|| for the case cr 2 > u\. In later sec- 
tions we consider materials with <j\ > 02 for which this 
bound is (see ref. |ll|), by interchanging the roles of the 
materials, 



<Tj = a 2 - 



l + (l + 2p)/3 12 -2( g Ci-p)/3i 2 2 
l + g /? 12 -(2 g Ci+p)/3? 2 : 



(8) 



where (3i2 — {cr\ — 02)l(p\ + 2o2). By way of mathemati- 
cal analogy these bounds also apply to to the effective di- 
electric, diffusion and magnetic permeability coefficients 
of composite materials. 



III. CORRELATION FUNCTIONS FOR THE 
GAUSSIAN RANDOM INTERFACE MODEL 



There is an extensive literature on the calculation of 
statistical correlation functions jl| . The case of the three 
point solid-solid correlation function has been consid- 
ered empirically p6|, p^|j4^ , [4^ ] , and theoretically for cellu- 
and spherical inclusions [|23|,pT 




lar materials 1 13 4J| and spherical inclusions |23| , |15| , |16| , [45[ 
to name a few. Here we take the interface between the 
phases to be defined by a level cut of a random field 
Now consider a Gaussian random field j/(r) 
(see section 0) and let the level sets y(r) = a 
define the interface (with the region y > a being phase 
1). Then the n point correlation function is given by the 
volume average, 

p„(ri,r 2 , . . . ,r„) =< H(yi - a) . ..H(y n - a) > (9) 

where H(y) is the Heavyside function and yi — y(jCi). 
p n is then the probability that the n points will lie in 
phase 1. For a macroscopically homogeneous isotropic 
material p n only depends on the distances ry = |rj — Tj\ 
between the points. Since volume and ensemble averages 



2 



are equivalent in such a medium pfl| we can use the lat- 
ter to evaluate eqn. (^). The joint probability density of 
Vi is, 



Pn{yW2 ■ --Vn) = 



, 1 eM-ly T G^y), (io) 

y/(2n) n \G\ 2 



where the elements of G are gij =< yiyj > j46|. The 
latter quantity we refer to as the field-field correlation 
function, 



k\Ti 



dk. 



(11) 



where p(fc) is the spectral density of the field. 

Berk |3^^^] and Teubner [B6| have derived the one 
point function (volume fraction), 



1 f°° 1 
p = -= \ exp(--t 2 )dt 

and the two point function, 



1 rsa 
P2{9ij) = ^ ex P 



dt 



P 



(12) 



(13) 



The three point function is calculated using the tech- 
niques described in ref. [p6| . The following identities are 
used g6), 

< exp(iy.w) >= exp(— -w T Gw) (14) 

and 

H(y-a) = — [ e - lw{ - y -^— (15) 

2lTl J c w 

where the contour C is directed along the real axis except 
near the origin where it crosses the imaginary axis in the 
upper half plane. This leads (after algebra) to, 



dpi _ dp 2 {gi2) 1 



= f exp(-~i 2 )<ft, 

7T JaF 12 1 



where, 



dgi2 

V 1 — 912 1 + 512 - 513 - 523 



(16) 



(17) 



VI + 5i2 v/|G| 

and \G\ = l-gf 2 - g 2 3 - 923 + 2 9i29i3923- Similar expres- 
sions can be derived for dp^ /dgi 3 and dpj /9g23- Defin- 
ing Aij = dpi /dgij we have therefore, 

pI (512, 5i3, 523) = 512 / Ai2(tgi2,tgi3, tg 23 )dt 
Jo 



+ 513/ A 13 (tg 12 ,tg 13 ,tg 23 )dt 



•f 

Jo 



+ 923 / A 23 (tgi 2 ,tg 13 ,tg 23 )dt. (18) 



The truncated three point correlation function pj is re- 
lated to the P3 by the expression, 

753(512,513,523) =pl (512,513,523) +PP2{gi2) 

+ PP2{g\3) + pp2{g23) - 2p 3 . (19) 

To examine the limit 7*12 — > (312 — > 1,523 ~^ 513) 
sct /(513) = P 3 (1,513,513) then df{gi 3 )/dg l3 = (1 - 
^P)dp2{gi3) I dgi3 and /(0) = implying p%(l, g 13l g 13 ) = 
(l — 2p)(p 2 (gi 3 )—p 2 ) as it should. (Similarly in the other 
limits) . The X-ray spectra of these materials can be cal- 
culated from p 2 1 22 |2^j3^] and hence they can be related 
to physical composites. Furthermore it has been shown 
that the surface to volume ratio is given by p^J36l], 



S_ 
V 



< k 2 >. 



(20) 



As the evaluation of the integrals in eqns. ( jl3| ) & (|l9| ) 
are computationally intensive it is useful to derive various 
approximations. Rigorous approximations for p 2 for the 
cases I a I -C 1 and |a| 3> 1 are derived in appendix [A|. A 
useful non-rigorous approximation to p\ 23 can be devel- 
oped by requiring that the approximation have similar 
properties to the actual function for r%j 3> 1 and sat- 
isfy the known consistency conditions in various limits 
||40f . Using the compact notation pj- = p 2 (g^ ) — p 2 and 

Pljk =P 3 (9i3,9ik,gjk) we have (r^ > 1), 



Pij x 9ij 

P\23 * 512513 + 512523 + 513523- 



(21) 
(22) 



Using this information, and including a higher order term 
for consistency {p 3 {r\ 2 , ri2, 0) = p 2 {r\ 2 )) we construct, 



T 
-Pl23 



1 - 2p 
2p(l-p) 
1 - 2p 



(Pl2Pl3 + P12P23 + P13P23 



2p 2 {\-pf 



T T T 
P\2P\3P23- 



(23) 



We note that this approximation has a maximum abso- 
lute error of O(10 -3 ) for the materials considered here. 
As such it is an order of magnitude better than a pre- 
viously suggested approximation |^7J p 3 (r \ 2l r 13, r 23 ) w 
P2(u)p2(v)/p where u and v are the smallest, and next to 
smallest, values of r 12 ,r 13 and r 23 . 



IV. DETERMINATION OF 

Actual calculations of the microstructural parameter 
Ci (eqn. (0)) have, to date, been for four classes of ma- 
terials. Cellular materials ]T^Jl^| , empirically measured 
physical composites [28], periodic arrays of spheres 
and materials with spherical inclusions. In the latter class 
the cases studied include: identical overlapping spheres 
]l5| , [il| (the IOS model), identical hard spheres Jigfl , and 
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poly-dispersed spheres |^5|,|| (many of these results are 
summarised in ref. 0). 

We now describe aspects of the computation of £1 for 
several spectra of the Gaussian random interface (GRI) 
model. It can be shown that Ci — 5 f° r P — 5 [ F^HF^l 
(see appendix |b|) and that Ci — 1 ~ C2 where £2 is the 
microstructure parameter associated with phase 2. As £1 
is dimensionless it must depend only on the ratios of the 
length scales associated with the spatial variables in the 
(dimensionless) correlation functions. That this should 
be so also follows from a simple dimensional analysis 
of the equations governing the electric field (no physical 
length scale is present). Henceforth we scale all spatial 
variables against a characteristic decay length without 
loss of generality. 

We consider three types of media generated from Gaus- 
sian random fields. The field-field correlation functions 
and their corresponding spectra are: 
Model I: 



9(r) 



vr 



p{k) = n- 2 ((1 - v 2 + k 2 ) 2 + Av 2 



(24) 
(25) 



Here v — 27r/ 1 /Z 2 with li the decay length of the field 
and I2 the characteristic domain size. When r » the 
correlation functions arisin g fr om this model are similar 
to those considered in refs. |22| (y = 0) and |^0| (y > 0). 
Note that this model has an infinite surface to volume 
ratio since < k 2 > diverges, however for computational 
realizations of the model the ratio is finite (see section 
0) and the model is well defined. However it is interest- 
ing to study (1 for this model to investigate the effect of 
interfacial roughness on effective properties. 
Model II: 



9{r) 
p(k) 



e 4 



(47T)i 



(26) 
(27) 



Model III: 



3 (sin fir — fir cos pr — sin r — r cos r) 
9( T ) = „ 3 /„ 3 i\ (28) 



p(k) 



47r(/i 3 - 1) 



r 3 (p 3 - 1) 
(H(p)-H(l)) 



(29) 



where p = fci/fco- Note that p(k) — > 8(h) as p — > 1 and 
the simple model used by Berk jBa ] is recovered. 

To perform the integration pwe use cylindrical co- 
ordinates (which damp the singularity at the origin), in- 
terchange the order of integration and exploit the r — s 
symmetry to give, 



I(w,cf>,u)dwd(f> P2(u)du (30) 




I(w,<p,u) 



pp 3 (r,s,t) -P2(r)p2(s) 
p 3 w sin 4> cos <fi 



(31) 



To elucidate the nature of the singularity in the integrand 
at w — 4> — we consider w, <j> (and hence r, s and t) 
to be small and assume the form g(r) ps 1 — ar n (where 
n = 1, 2 in accord with models I-III). Now for p —\ the 
numerator of the expression for / is given by, 



''{git)) sin 1 (g(r)) sin 1 (g( S )) 
8tt 4tt 2 



y/2aw n cj) n , (32) 



n 



- sin 2(f>u. Thcrc- 
i and numerical 
i An 



m- 



where we have used the results sin 1 (g(r)) w ? — \/2ar 
r = w cos <j>, s = wsxncf) and t = w^/T 
fore I(w,(f>,u) ~ (uxf))^- 1 with p = 
analysis shows this scaling also holds for p =/= 
tegration rule which takes the singularity into account 
is employed. Note that for finite surface area to volume 
ratios n = 2. The number of abscissae in each of three 
integration ranges was increased until the third signif- 
icant figure in the estimation of (1 remained constant. 
The integration method was tested on the known corre- 
lation functions for the IOS model [ p3|Jl5| . The results 
are in exact agreement (to the reported 3rd significant 
figure) with those of Torquato and Stcll Jl],[l5| and agree 
to the second significant figure with the results of Berry- 
man pl|. The calculation of the correlation functions 
in the integrand is done using a combination of iterative 
quadrature rules |50| , the asymptotic results presented in 
appendix |X| and the non-rigorous approximation for the 
truncated three point function (p3|). The latter is used 
whenever 1 — g < 1CP 5 since in this case the functions Fij 
exhibit large derivatives and the quadrature rules con- 
verge too slowly. The accuracy sought in the application 
of each approximation is O(10 -6 ). 

TABLE I. The microstructure parameter fi for various 
GRI models. Prior results for the identical overlapping 
sphere (IOS) model |l|[l5| and the symmetric spherical cell 
(SSC) model JL3IJ7I are included for purposes of comparison. 



Prior Results 



GRI model 



SSC IOS 



1^ = I v = 10 



II 



III (j, = 1.5 



0.01 
0.05 
0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
0.95 
0.99 



0.01 
0.05 
0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
0.95 
0.99 



0.056 
0.114 
0.171 
0.230 
0.290 
0.351 
0.415 
0.483 
0.558 
0.604 
0.658 



0.269 
0.293 
0.319 
0.366 
0.411 
0.456 
0.500 
0.544 
0.589 
0.634 
0.681 
0.707 
0.731 



0.193 
0.217 
0.248 
0.309 
0.372 
0.436 
0.500 
0.564 
0.628 
0.691 
0.752 
0.783 
0.807 



0.099 
0.160 
0.210 
0.291 
0.363 
0.432 
0.500 
0.568 
0.636 
0.709 
0.790 
0.840 
0.901 



0.060 
0.105 
0.150 
0.237 
0.324 
0.411 
0.500 
0.588 
0.675 
0.763 
0.850 
0.895 
0.940 



with 
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where, 
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0.0 [ 



0.0 



0.2 



0.! 



1.0 



0.4 0.6 
P 

FIG. 1. The micro-structure parameter vs. p for the 
GRI models (equations @-@) and the IOS model @jl5). 

The results for each of the models is presented in table 
| and plotted in figure [j] along side the results for the IOS 
model. Several comments on the qualitative relationship 
between £i for the GRI model and prior calculations can 
be made. Consider the expansion, 



Ci 



E 

i=0 



&iP'. 



(33) 



For a general class of materials with spherical inclusions it 
has been shown that eo = . If the inverse of these 
materials is considered (so that the correlation functions 
refer to the material surrounding the inclusion) or el- 
lipsoidal inclusions are considered Jl7|] then eo > (see 
appendix |^). For the case of symmetric cell materials 
eo = M £ [0, 1] where M = for spherical cells and 
M = 1 for plate like cells ]7],[l3| ■ Another interesting fea- 
ture of Ci is that it is observed to be linear with p over 
a wide range Indeed for the symmetric cell model 

ei = (1 — 2AI) and e2, e 3 . . . = 0. By inspection of hgure 
[l] we see that eo > in qualitative agreement with the 
results for non-spherical inclusions which will occur in 
the GRI models. Also note that £i is very similar to the 
results for the symmetric cell model for some < M < 1 
over a wide range of p. This discussion demonstrates 
the success of simple models of random media to capture 
the qualitative features of £i for the amorphous materials 
considered here. 

Calculations of the related microstructure parameter 
r/i which arises in bounds on the elastic properties of 
random composites are reported in appendix M. 



V. GENERATION OF FIELDS 

For computational purposes we consider a T-periodic 
Gaussian random held |16j with a maximum wavenumbcr 
K = 2nN/T, 



N 



N 

E 

--N m=-N n=-N 



N 

E 



(34) 



2n n- • ,n 

iron = fri 11 + TO J + nk )- 



(35) 



For y K real we require ci :m , n = c_;,_ m ,_„ and as < y K >= 
we set co,o,o = 0. For reasons which become clear be- 
low we also take c/ m „ = for ki mn = |k; rmi | > K. With 
cimn = aimn + ibimn, CLlmn and bi mn are random indepen- 
dent variables (subject to the conditions on c/„ m ) with 
Gaussian distributions such that < ai mn >—< bi mn >= 
and 



=< bf„ 



>= 



1 



Pk (kimn) 



(36) 



with p K (k) the spectral density. The field-held correla- 
tion function g K is given by 



(37) 



c(ri 2 ) = < 2/K(ri)y4r 2 ) > 

N N N 

= E E E< c imnc^ > e lki — ( ri ~ r2 > (38) 



—N —N -N 

\^ PK (k) 8h l^-^ dk. 

fc|ri-r 2 | 



(39) 



The last integral is obtained by taking TV and T large, us- 
ing equation ( |36"|) and recognising that the summation is 
the approximation of a triple integral. Following a trans- 
formation to spherical co-ordinates we integrate over the 
angular variables to obtain j39|). Since g K (0) = 1 we 
define, 



r K 

p K (k) = p- 1 p(k), P= 47Tfc 
Jo 



' p(k)dk 



(40) 



where the p(k) arc defined in section [II. If in addition 
we take K — > oo then the conventional correlation func- 
tion (and spectral density) are recovered. The Fourier 
expansion (|34| ) is evaluated using a FFT. 

Consider materials derived from the field y K (eqn. (|3~i|)) 
as discussed in section [II. Cross sections of the media for 
four different variants of the models are plotted in figure |^ 
(a)-(d). The large scale structure of the interface is deter- 
mined by the terms in the expansion with small k while 
the small scale structure (ripples on the surface) are de- 
termined by the terms where k is large. A physical mate- 
rial will naturally contain a finite cutoff wavelength cither 
imposed by the molecular size or by the manifestation 
of surface tension at the phase boundaries. This wave- 
length will then dictate the grid resolution necessary to 
properly resolve the structure for simulational purposes. 
Conversely for this study the discretisation Ax = T/M, 
where M is the number of grid points in the x direction, 
will restrict the choice of K. Thus Ax <C \ m in = 2tt/K 
or N M. Another constraint 
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(a) I: v=0, K=18, T=4ji (b) I: v=10, K=18, T=4jc (c) II: K=18, T=4ji 



(d) III: u=1.5, T=4ji 




(e) I: V=0, K=8, T=4ji 



(f) I: v=10, K=32, T=3i 



(g) II: K=8, T=47t (h) III: (1=1.5, T=8n 

FIG. 2. Cross sections of the models generated at the same scale (a)-(d) and at the scale to be used in 
the simulations (e)-(h). The volume fraction p — 0.5 and each of the fields are generated using the same 
random number seed to clearly show the differences amongst the models. The parameters v, fi, T and K 
are discussed in the text. 



is that T 3> Id where Id is an effective decay length of 
the field defined by g(ld) — e _1 . This must be so to 
ensure that the edges of the sample are uncorrelated to 
properly simulate an infinite random medium. For the 
approximation involved in determining g K (eqn. (^)) to 
be accurate would require T 3> 27r, however in this study 
this constraint was found to be of less importance and is 
ignored. 

The computational parameters used for each of the 
four variants of the models are given in table [n]. Cross- 
sections of each of the materials are plotted in figures || 
(e)-(h) and the interface for models II and III are plotted 
in figures || and ||. To properly account for the effect of 
K for the computational models d must be recalculated 
using g K (r)- This is done using a look up table gener- 
ated by numerical integration of ( fj^ ) and an asymptotic 
expansion for large r (where the integral is costly to eval- 
uate), 



9k{t) = -pg{r) 



AiTKp K (K) cos K r-~- 



O I ^3 



(41) 



The function g K {r) for model I is plotted along with 
direct measurements of < y(ri)y(r 2 ) > in figure g. The 
agreement is seen to very good. The values of £i for the 
computational models are presented in table III. The 
effect of K on £i can be seen most clearly for Model I 
where the surface area to volume ratio becomes arbitrar- 
ily large as K increases. The interfacial smoothing effect 
of imposing a finite cut-off (K) is shown in figures ^ (a) 
(K = 18) and (e) (K = 8). 



TABLE II. The parameters used to gener ate the compu- 
tational materials discussed in sections 
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Model 


T 


N 


K 


I 


v = 


4-7T 


16 


8 


I 


v = 10 


71 


16 


32 


II 




4tt 


16 


8 


III 


H = 1.5 


8tt 


6 


1.5 




FIG. 3. A plot of the interface y K {v) = (p = 0.5) 
for Model II. The parameters used to generate the field are 
K = 8, T = 47T. Note that the large scale structure of this 
model is similar to that of model I (y — 0) as can be seen by 
comparing figs. || (e) & (g). 
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TABLE III. The microstructure parameter £1 for the com- 
putational materials listed in table |l[ The parameter K is 



FIG. 4. A plot of the interface y K (r) = (p = 0.5) 
for Model III. The parameters used to generate the field are 
(j, = 1.5, T — 8ty. 

This has the effect of decreasing (increasing) the mag- 
nitude of for p < | {p > |) (compare tables | and |l|). 
For small (and large) volume fractions this can be quali- 
tatively explained by the fact that the inclusion phase of 
the model will be much rougher (and hence less 'spher- 
ical') for K = oo than for K = 8 (See appendix Q). £i 
for model II is unchanged for finite K to the accuracy 
calculated here. 

VI. SOLUTION OF THE LAPLACE EQUATION 

There are a variety of different methods of simulating 
the effective conductivity of an inhomogeneous medium. 
These include direct solution of the partial differential 
equations governing the potential <fi [ pl| , Brownian mo- 
tion algorithms P-|i"T|, Fourier methods |32j and other 
techniques |l2[. 




o.oo 



0.10 



0.20 



r/T 



0.30 



0.40 



0.50 



FIG. 5. The computational field-field correlation func- 
tion gx(r) (eqn. ([ji])) for model I (y = 0, K = 8). The exact 
form (ie. K = oo) given by eqn. (^i|) and the asymptotic 
form, eqn. ([fl|), are included for purposes of comparison. 
The data was measured from a single sample. 



the maximum 


wave number in 


the Fourier Series. 


The rela- 


tion = 1 


- £i(l - p) can 


be used to determine (i for 


P> 1/2. 








P 


lu = 


I v = 10 


II 




K = 8 


K = 32 


K = 8 


0.01 


0.155 


0.085 


0.099 


0.05 


0.214 


0.136 


0.160 


0.1 


0.258 


0.182 


0.210 


0.2 


0.326 


0.265 


0.291 


0.3 


0.387 


0.344 


0.364 


0.4 


0.444 


0.422 


0.432 



We solve Laplace's equation with the charge conser- 
vation boundary conditions discussed in section || in a 
cube of side length T using M 3 nodes. A potential of (pi 
and 4>o are applied on the faces z = and z = T and 
periodic boundary conditions are imposed on the four 
faces parallel to the direction of the current to model a 
sample of infinite extent in the x and y directions. A 
finite difference scheme is used to approximate the field 
and is solved using a conjugate gradient method. In ap- 
pendix [S] we discuss the efficient implementation of the 
algorithm on a parallel computer. The z components of 
the current and the field are then used in equation (^) 
to determine (cr e )M- The convergence criterion of the 
CG solver is decreased until the third significant figure 
of (cr e )i\i remained constant. To estimate the continuum 
value of tr e we assume that a e « (cr e )M + a\M~ x and fit 
a line (using least squares) to several values of (a e )M vs. 
M _1 . The intercept of this line with the axis M~ l = 
then provides a e - 

Before proceeding to the random media we simulate 
the effective conductivity of a regular array of spheres of 
conductivity ct\ = 10 in a matrix of conductivity a% = 1. 
Exact results for this model have been calculated by 



McKenzie et al. 
previous authors 



and the model has been used by 
10,|12| to test the accuracy of their al- 
gorithms. For computational purposes the array contains 
four spheres in the z direction (using six spheres changes 
a e by less than 1%). The values of (<J e )M for increasing 
concentration are plotted along with the lines of best fit 
used to estimate a e in figure ^|. The graph demonstrates 
the necessity of extrapolating the data to M _1 — > 0. The 
results for er e are presented in table IV. The error is less 
than 1% for p < 0.4 but increases to around 3% at p = 0.5 
near the percolation threshold p c « 0.52. For the random 
media it was found that computing (a e )M at more than 
two values of M did not significantly alter the estimation 
of a e . 

For the random media we must consider how to assign 
the conductivity of a bond lying between two nodes which 
lie in different phases. Let yij and a^j be the respective 
values of the field and the conductivity at two 
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p=0.5 




0.030 



0.000 0.010 0.020 

1/M 

FIG. 6. The calculated values of (ce)*/ (solid circles) for 
varying volume fraction (p) and the exact data (open circles) 
|33| for the regular array of spheres. M is the discretisation 
used in the finite difference scheme. The lines are linear 
least squares fits of the data over the range of M _1 where 
the data is approximately linear (generally M _1 < 0.016). 
<7 C is then then given by the intercept of the lines with the 
axis M' 1 = 0. 

such neighbouring nodes. There are three obvious ways 
of determining cr^ . Defining a = (a — yi)/ (yj — yi) we can 
choose Gij = acrj + (1 — a)<Ji (as if the portions of the vol- 
ume element associated with the bond are like conductors 
in parallel), — (a/<Tj + (1 — a)/o"i) _1 (as if in series) 
or 0V, = o i or er 2 as (y t + yj)/2 > a or (y t + yj)/2 < a 
(a simple field average). In figure we show the effect 
of using these rules for two samples of material I gener- 
ated with N = 4, 16. For a given discretisation (M) a 
large difference in a e occurs depending on the rule em- 
ployed. However extrapolation to M _1 = demonstrates 
remarkably well that the choice is immaterial. As the 
simple averaging rule provides the least error for given 
M it will be used. Finally, for a given volume fraction 
we use the bisection method to calculate the value of 
the level cut parameter a' such that < a > = po\ + qo^ 
(where <> refers to bond averaging). This substantially 
reduces the statistical fluctuations in a e compared to us- 
ing^ the theoretical value of a determined from equation 

TABLE IV. Comparison of the extrapolated finite differ- 
ence simulations with the exact results |33| for a simple cubic 
array of spheres of conductivity o\ = 10 in a matrix of con- 
ductivity o"2 = 1. The Brownian motion simulations of Kim 
& Torquato jjl| and the simulations of Bonnecaze & Brady 
are also included. The spheres touch at p = 7r/6 ~ 0.52. 





Exact 


Finite Diff. 


Relative 


Kim- 


Bonnecaze- 


p 


Results 


This work 


Error 


Torquato 


Brady 


0.1 


1.24 


1.25 


0.8% 


1.24 


1.24 


0.2 


1.53 


1.52 


0.7% 


1.53 


1.53 


0.3 


1.89 


1.90 


0.5% 


1.89 


1.87 


0.4 


2.36 


2.37 


0.4% 


2.36 


2.29 


0.5 


3.11 


3.19 


2.6% 


3.13 


2.80 
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FIG. 7. The values of (cr e )M for sample Gaussian media 
using the three different rules for assigning a conductivity to 
each bond (see text). Note that large variations in (a e )M oc- 
cur for given M if different rules are employed. However the 
choice is seen to be immaterial when the data is extrapolated 
to M" 1 -> 0. 

VII. SIMULATION RESULTS 

We simulate the effective conductivity for the four GRI 
models listed in table |H| for a range compositions. As we 
are dealing with finite sized samples we report cr e as the 
average over a number of different realizations of the ma- 
terials. Error bars, which represent 95% confidence limits 
on the results, are equal to twice the standard error. The 
samples are examined at three different contrast values; 
01,2 = 10, 1, gjL,2 = 50,1 and 01,2 = 1,0. Previous au- 
thors Q; 
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5J| have considered media with <7i,2 = 00, 1 
but this is not possible using the methods discussed here. 

The results for the case 01,2 = 10, 1 are presented in 
table |v| and plotted in figure | for p £ [0.1,0.9]. To 
obtain the data five samples of each model with discreti- 
sations using M — 64 and M = 96 were considered. The 
results show little variation in a e (for fixed p) amongst 
the four different media. The maximum relative differ- 
ence of 4.2% occurs at p — 0.6 between a e for model I 
(v = 0, K = 8) and model III (fi = 1.5). As the dif- 
ferences are relatively small we restrict further attention 
to the latter two materials. The bounds calculated from 
equations (§) and j]8|) using £1 from tables Q and |l| are 
and plotted along with the simula- 
, for p e [0.2, 0.8], and figure [l0| for 
p e [0.86,0.96]. The latter figure illustrates very clearly 
that the upper bound discriminates between model I and 
III. 

For the case <j\p = 50, 1 the results for model I {y = 0, 
K = 8) and model III (/i = 1.5) are reported in ta- 



presented in table K 
tion data in figures 



ble VII and plotted along with relevant bounds in figure 
[n]. Again 5 samples of the media were considered with 
M = 64 and M = 96. Figure [fl] shows very pronounced 
differences between the IOS model and GRI mo dels. The 
results for the case <Ti 2 = 1,0 are given in table VIII and 
plotted along with the upper bound (the lower bound 
vanishes) in figure [l2| Five samples at discretisations 
M = 48 and M = 64 were considered. 
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TABLE V. Simulations of a e for different GRI models 
with 0-1,2 = 10, 1. The data for the IOS model @ is in- 
cluded for purposes of comparison. The notation x(y) im- 
plies cr e = x ± y.lCT 2 with the error bars defined as 95% 
confidence limits a . 
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4.21(4) 
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5.12(2) 
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5.24(6) 
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6.23(1) 
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6.33(6) 
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7 41 (1 ") 


7.39(1) 
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8.13(3) 






8.20(4) 




0.88 


8.39(3) 






8.46(4) 




0.90 


8.65(2) 
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9.45(1) 






9.48(1) 
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o III: (1=1.5 
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FIG. 8. The data for four different GRI models with 
conductivities <ti = 10 and o\ — 1. The data is normal- 
ized against the BM upper bound (eqn. (g)) with (i = p to 
highlight the differences between the models. The curves are 
spline fits of the data. 



TABLE VI. An example of the bounds calculated using 
equations (^) & (^) and the data in tables 
(Model III). The conductivity contrast is a\,% 



(Model I) & | 
10,1 





I K=8, v 


= 


III fi = 


1.5 


p 


<Ji 






cr u 


0.1 


1.290 


1.437 


1.268 


1.372 


0.2 


1.660 


1.994 


1.618 


1.905 


0.3 


2.123 


2.654 


2.073 


2.576 


0.4 


2.702 


3.414 


2.663 


3.371 


0.5 


3.423 


4.273 


3.423 


4.273 


0.6 


4.316 


5.229 


4.385 


5.269 


0.7 


5.414 


6.283 


5.570 


6.348 


0.8 


6.742 


7.434 


6.961 


7.501 


0.9 


8.302 


8.677 


8.490 


8.720 
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FIG. 9. Simulations and bounds for the effective conduc- 
tivity of two models of random composite media generated 
from the GRI model for the case ai^ = 10, 1. The circu- 
lar symbols represent data for the IOS model calculated by 
Kim & Torquato |nj and the bounds for the IOS model were 
evaluated by Torquato & Stell Jl5[ . 




FIG. 10. Simulations and bounds for the effective con- 
ductivity of two models of random composite media gen- 
erated from the GRI model for the case 1X1,2 = 10, 1. The 
upper bound discriminates between the models in this range 
of p. 



0.90 0.92 0.94 0.96 








FIG. 11. As in figure ^except o\ — 50 and a% = 1. 

TABLE VII. Simulations of cr e for the different Gaussian 
random media with a\ = 50 and C2 = 1. The IOS model 
data is from ref. fill. The error bars define 95% confidence 



limits. 

I^ = 0,_ft: = 8 HI /j = 1.5 IOS 

(U 1.88±0.08 1.69±0.06 

0.2 3.68±0.2 3.45±0.2 2.16 

0.4 10.4±0.5 11.1±0.3 6.44 

0.6 21.0±0.4 22.5±0.3 15.2 

0.8 34.4±0.4 35.6±0.2 30.7 




FIG. 12. As in figure |^ except a± = 1 and 02 = 0. 



TABLE VIII. Simulation of a e for different Gaussian ran- 
dom media with o\ = 1 and 02 = 0. The IOS model data 
is from ref. jnj. The notation x(y) implies a e — x ± y.l0~ 3 
with the error bars defined as 95% confidence limits a . 
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I, v = 0,K = 8 


Ill /j = 1.5 


IOS 


0.1 


0.002(2) 


0.000(0) 


0.022 


0.2 


0.027(6) 


0.025(2) 


0.076 
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0.080(3) 


0.092(6) 


0.160 


0.4 


0.166(12) 


0.185(10) 


0.248 


0.5 


0.264(10) 


0.294(12) 


0.346 


0.6 


0.383(9) 


0.420(11) 


0.461 


0.7 


0.521(12) 


0.553(6) 


0.593 


0.8 


0.670(7) 


0.697(7) 


0.714 


0.9 


0.830(3) 


0.847(2) 


0.855 


a y = 


implies y < 0.5. 







There are several qualitative trends in the data which 
can be commented on. Note that for first two contrasts 
considered a e of the GRI models is greater than that of 
the IOS model over the entire range of p. At low volume 
fractions this can be attributed to the fact that the inclu- 
sions of the GRI models are qualitatively less spherical 
than those of the IOS model (see appendix |c|) . This can 
be clearly seen in figures [13] and [l4| where the inclusions 
are plotted for each of the GRI models at p = 0.07 (the 
IOS model will containing predominantly spherical inclu- 
sions at this volume fraction). At high volume fractions 
the situation is reversed; the matrix phase of the IOS 
model is extremely ramified and hence a e is lower. Simi- 
larly near p = 0, 1 the small differences in a e for Models I 



and III can be explained by the fact that the latter model 
has more spherical inclusions (compare figures |l^ & [TI]) . 
This behavior is consistent with the relative variations in 
Ci for each of the three models as discussed in appendix 
|Cj. For mid-range p the differences between the IOS and 
GRI models correspond to the fact that the more highly 
conducting regions of the latter are generally better con- 
nected than those of the IOS model. Again this difference 
can be anticipated from the relative behavior of the pa- 
rameter £i for the two classes of models. However this 
is not necessarily always so as can be seen by comparing 
the respective values of Ci (tables | & III) and <j e (table 
0) for Models I and III at p = 0.4. In this case a\ < a" 
but Cl > C? 1 - 



in 
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FIG. 13. A plot of the interface y K (r) = 1.48 (p = 0.07) 
for Model I. The parameters used to generate the field are 
v = 0, K = 8, T = 4n. 

Note that the simulation data for the IOS model in the 
case <7i,2 = 1,0 were obtained for insulating spheres in 
a matrix of unit conductivity. Therefore the microstruc- 
ture of the conducting phase is generally better connected 
than the GRI models and the arguments pertaining to 
the qualitative variations in a e used above are reversed. 

In all of the above cases the data lies between the ap- 
propriate bounds. Furthermore the upper bound is seen 
to provide a good estimate for a e for the GRI models 
for p > 0.7. This has been observed previously ]i], p5| , pl| 
for both high and low p (the latter case has not been 
considered in detail here). This fact provides evidence 
of the near optimality of bounds in these regions and a 
good confirmation of the techniques used in this paper 
(see figure [lO] especially) . 

We have calculated the percolation threshold for Model 
I and III using the algorithm of Skal et al. ||(J. In con- 
trast to previous results we found p c to exhibit finite 
size effects and to depend on the spectrum. Averaging 
p c over 10 fields at M = 20,32,48,64 and extrapolating 
to M = oo we found p c as 0.07 (a c ~ 1.47) for Model I 
and p c » 0.13 for Model III (a c ~ 1.13). In the absence 
of theoretical percolation results for the GRI model it is 
interesting to discuss the threshold in terms of the tran- 
sition of the structures from elliptic to hyperbolic as p 
increases. The average Gaussian curvature [Kq) for the 
GRI models is given by |3(| , 

< K G >= ~ < k 2 > (a 2 - 1). (42) 

Therefore the nature of the interface is dominated by in- 
clusions of positive curvature (eg. ellipsoids) for \a\ > 1 
and is predominantly hyperbolic (eg. bicontinuos) for 
\a\ < 1. The fact that a c > 1 indicates that connecting 
structures persist below the elliptic/hyperbolic transition 




FIG. 14. A plot of the interface y K (r) = 1.48 (p = 0.07) 
for Model III. The parameters used to generate the field are 
p = 1.5, T = 8tt. 

as would be expected since < Kg > is an average quan- 
tity. 

Finally it is interesting to discuss the surprisingly small 
variation in a e (and (V) amongst the GRI materials. As 
can be seen in figures || (a)-(d) these materials appear to 
be very different when viewed at the same scale. How- 
ever the major qualitative differences arc related to the 
effective decay length, and when the materials are viewed 
at a scale proportional to this length they appear to be 
remarkably similar as shown in figures ^ (e)-(h). Equiva- 
lently if the length parameters are retained in eqns. (pif)- 
( p9| ) they can be tuned to achieve the latter group of 
figures at the same scale (without effecting er e ). The re- 
maining smaller qualitative differences account for the 
variation in a e observed. 



VIII. CONCLUSION 

We have investigated the effective conductivity of a 
two phase random composite material using bounding 
techniques and direct simulation. Our calculations of d 
increase the classes of composites to which the bounds 
can be applied, while the simulation data can be used 
to assess predictive theories ||6[ and be compared with 
higher order bounds. The results also pertain to a variety 
of other effective properties of amorphous composites as 
discussed in section [n} 

The bounds encompass all of the simulational data and 
the upper bound yields a reasonable estimate of er e for 
p > 0.7 (ci > (T2). Reasonably large differences in <j e 
and Q\ are observed between the amorphous GRI mod- 
els and the IOS model. This highlights the importance 
of incorporating microstructure effects in the calculation 
of the effective properties of composite materials. Con- 
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versely there is relatively little variation in a e amongst 
the GRI models as the major qualitative microstructural 
differences between the models are related to an effective 
decay length upon which a e is necessarily independent. 
We expect that other properties of such composites where 
no intrinsic field length scale is present (eg. the elastic 
bulk and shear moduli) will show similar behavior. 

It is clear that the Gaussian random interface model 
discussed here can serve as a useful 'model' amorphous 
medium in the study of the effective properties of ran- 
dom composites. Furthermore the bounds and simula- 
tions can be related to physical composites by experi- 
mentally relating such systems to one of the theoreti- 
cally known models. This could be done by comparing 
the spectra obtained from small angle scattering studies 
with that obtained from the 2-point correlation functions 
of each model. Or, more simply, by comparing images of 
the models with electron micrographs. Although such 
schemes are only approximate, our results indicate that 
the fine microstructural details are relatively unimpor- 
tant in both the calculation of £1 and the simulation of 
a e . 

We note that the GRI model can be extended to the 
case of membranes and foams f3^]56|] and that higher 
order correlation functions can be calculated for use in 
more precise bounds. Random walker algorithms can be 
utilised to investigate the often studied scaling properties 
of a e near the percolation threshold. 



For the case a>0we have (using successive integra- 
tion by parts), 



2ira 



1 + 



(-l)"1.3...(2n 



n=l 



(A3) 



Special care must be taken to determine a practical ex- 
pansion for p2- A simple application of Watson's Lemma 
p7[ yields a solution which does not possess the correct 
limiting behavior as g — > and gives two different expan- 
sions for the cases g — 1 and g < 1. The first problem 
is dealt with by appropriately partitioning the integral, 
while the second necessitates a further transformation of 
variable, followed by an expansion in scaled parabolic 
cylinder functions (see for example ref. J58|). Thus we 
write 



P2i9) = h (/_r/J cxp (-it* 



dt 



+ tj y/T~P 



■P ■ 

(A4) 



In the second of these integrals we make the substitution 
v = 1/(1 + 1) - 1 to give 



-dv 



2^TT Jo (V + 1)(V + |)2 

which can be expanded using Watson's Lemma, 



(A5) 
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APPENDIX A: ASYMPTOTIC FORMS OF P 2 



The two point function (eqn. (|13|)) can be expanded in 
powers of a to yield, 



1 



P2 ( 5) = --arcsin ff +^— £ 

n=l 



with, 



" (-l)"a 2 " q» 
2tt ^ 2 n n\ 



P 2 (Al) 



2 / fl-g 
a n = II' 



2n- 1 



1+5 



(A2) 



and do = 0. This expansion converges rapidly for small 

Id. 



2tt 



1 2 7/1 



(A6) 



This is just the exp ansion of p 2 which can be cancelled 
from equation (A4). The remaining integral is put in a 
standard form by the substitution v = 1 /(!+£) — l/(l+g), 



P2(g) 



2** Jo (v + 0-)(v + $te)i 



dv. (A7) 



Note that the nature of the singularity of the integrand 
changes order as g — > 1 (this makes it impossible to gen- 
erate an expansion for the full range of g using Watson's 
Lemma). To develop a uniform expansion for large a 
valid near g = 1 we make the further substitution, 



1 



igU, 



$ 9 = 



1-ff 
1+3' 



to give 



2tt 



-a 2 (±u 2 +5 g u) 







6gU 



1 

1+9 



du. 



(A8) 



(A9) 



In the usual way the non-exponential component of the 
integrand can be expanded in powers of u and integrated 
term by term to give, 
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P2(g) R 

-f 

where, 
T„(z) 



1 ( 2 ^ 

(l-2g)(l+ g ) s 
2a 3 



^3 ( V) 



(A10) 



- x ds = T(n)ei z2 \J(7 



(All) 



U(a, 2) is a parabolic cylinder function p9| . Two simple 
checks can be made on this expansion. For g — 1 w e hav e 
T n (0) = 2"/ 2 ~ 1 r(§ ) which, when substituted in (|aTq|) , 
gives the expansion of p. For g < 1 we again employ Wat- 
son's Lemma to determine the asymptotic expansion, 



(-l)'T(n + 2i) 



n+2j 



(A12) 



Now taking g = in ( A10) and using this expansion gives 
the asymptotic form of p 2 as it should. The expansions 
for the case a -C (p sw 1) are simply, p = 1 — AEi 
and P2Q7) = 2p — 1 + AE 2 where AEi and AEg are the 
asymptotic expansions given by @ and ( |aTc|) respec- 
tively with a replaced by \a\. Note that Berk |35| , |37t has 
derived a formal series representation of P2 valid for all 
a, however the convergence of the series is slow for g w 1 
and not guaranteed at g = 1. 



APPENDIX B: PROOF OF Ci = \ FOR P = §. 

In Brown's formulation the parameter £1 arises as the 
limit as e — > of the integral, 



Ci 



2pq 



°dr 



duP 2 (u)f{r,s,t) (Bl) 



where / is the term in brackets of eqn. d?]). Now by taking 
p = \ in eqn. © we have / = j^(t)/2 - 2^(r)^(s). 
Note that the integral of the second term vanishes since 
it does not depend on u and J, P^iujdu = 0. Therefore 



after making the substitution t 



2sru eqn. 



(Bl) becomes, 



Ci = 9 / dr ds 



r+s 



pl(t)h(r,s,t)dt 



(B2) 



where h(r, s, t) = t(srT 4 (f(i 2 



2\2 ^2+2 



t 2 ). Inter- 



changing the order of integration results in, 



f2e />oo 

u = !»i / + y 2c ) vi\ 



rt+s 
ds I hdr 



9 / Pi (t) 

he 



t/2 ft-s 

1 ' 1 1 ' ds hdr, 



(B3) 



and carrying out the straight forward integrations over r 
and s leads to, 



Cl = 7o pl > {t) \ 2,' 24," 



dt. 



(B4) 



Now by taking e — > and using the fact that p^(O) = 4 
a _ 1 

= 1 - 2- 



gives Ci - 1 



APPENDIX C: RELATIONSHIP BETWEEN a E , <i 
AND SHAPE AT LOW P 

Consider the small concentration approximation J6C| ] 
to a e for the case of randomly distributed and oriented 
spheroids with axial ratios A, A, 1 — 2A, 

<r e = V2 + -p((Ti -a 2 )z(a u <r 2 ,A) + 0{p 2 ) (Cl) 



where, 



2ct 2 



(T 2 



2 + A(<7i - (7 2 ) £7 2 + (1 - 2A)(CTi - £7 2 ) 



(C2) 



With A S [0, 1/2] it can be easily shown that z has a 
unique minimum at A = 1/3 (spherical inclusions) and 
is monotonically increasing as \A— 1/3| increases. There- 
fore with a 1 > £72, £7 e will be higher the lower the spheric- 
ity of the inclusions (and conversely for cr 2 > £?i). The 
same argument should qualitatively hold for arbitrary 
shapes. To see how this relates to Ci we match the terms 
of the expansions of eqn. (|Cl|) and eqn. (@) to order p and 
(cti-ct 2 ) 3 which gives flJlJCi = (l-3A) 2 + 0{p). Thus 
Ci will be higher for less spherically shaped inclusions. 



APPENDIX D: IMPLEMENTATION OF THE 
FINITE DIFFERENCE SCHEME 

The finite difference scheme (see for example ref. | p(j| ) 
for the equations and boundary conditions discussed in 
section || leads to a system of simultaneous equations for 
the value of the potential at each of the interior nodes 
(including those on lateral faces if we define 4> to be pe- 
riodic in the x and y directions). For each such node u 
we have 



.0 = 



(Dl) 



where nn is the set of nearest neigbours of node u and 
a uv is the conductivity of the bond lying between nodes 
u and v. Conventionally these equations are cast as a 
matrix equation with <p and the boundary conditions as 
1 dimensional matrices (vectors) 1 50 . On a parallel com- 



puter it saves coding and implementation time to retain 
the potential 4>i,j,k as a 3D matrix. Define A as the op- 
erator which performs the operations defined on the left 
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hand side of (Dl) for interior nodes and A<p = <j> on the 
nodes where Dirichlet conditions are to be applied. Also 
define b as a 3D matrix containing the boundary condi- 
tions on the field (bij : i = 4>i, = <fio for ah hj an d 
6 = elsewhere). Then solving the system of equations 
for is equivalent to minimising || A<p—b \\2, which can be 
done using a conjugate gradient method which handles 
vectors of general dimension. 



APPENDIX E: CALCULATION OF THE 
MICROSTRUCTURE PARAMETER jji 

Three point bounds have also been derived for the elas- 
tic bulk and shear moduli [|)-0|6^] which can be expressed 
in terms of £1 and an additional parameter, 

5 150 f°°dr f°°ds f 1 , , , 
?/1 " Cl + Wo - Jo tL^M* 



21 

P3(r,s,t) 



r Jo - s 

P2{r)p2(s 
P 



(El) 



We have calculated 771 for several different GRI models 
and tabulated (see table IX) it along with data for the 



IOS model. Qualitatively the results are similar to those 
for Ci and we expect the differences in the effective shear 
and bulk moduli amongst the GRI models to be similar 
to those observed for the conductivity case. 

TABLE IX. The microstructure parameter rji which 
arises in bounds on the elastic bulk and shear moduli (see 
appendix JE ) for selected Gaussian media. Data for the IOS 
model 0,q2 | is included for purposes of comparison. 



p 


IOS 


Si 


I v = 

K = 8 


III (i = 1.5 


0.1 


0.075 


0.276 


0.213 


0.106 


0.2 


0.149 


0.333 


0.291 


0.197 


0.3 


0.224 


0.388 


0.362 


0.294 


0.4 


0.295 


0.444 


0.432 


0.396 


0.5 


0.367 


0.500 


0.500 


0.500 


0.6 


0.439 


0.556 


0.568 


0.604 


0.7 


0.512 


0.612 


0.637 


0.706 


0.8 


0.583 


0.667 


0.709 


0.803 


0.9 


0.658 


0.724 


0.787 


0.894 
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